#include <GL/glew.h>
#include <GL/freeglut.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include "Volume.h"
#include "FileManager.h"
#include "png.h"

typedef unsigned int uint;
typedef unsigned char uchar;

#define MAX_EPSILON_ERROR 5.00f
#define THRESHOLD         0.30f

uint width = 512, height = 512;

float * viewRotation		= NULL;
float * viewTranslation 	= NULL;
float 	invViewMatrix[12];
float 	density 		= 0.05f;
float 	brightness 		= 1.0f;
float 	transferOffset 		= 0.0f;
float 	transferScale 		= 1.0f;
int 	linearFiltering 	= 1;

GLuint 	pbo 			= 0;     // OpenGL pixel buffer object
GLuint 	tex 			= 0;     // OpenGL texture object

unsigned int timer = 0;

volume_t volume;
uint * 		d_output;

#define MAX(a,b) ((a > b) ? a : b)

/*
 *------------------------------------------------------------------------------------------------------------------------------------------------------------------
 *------------------------------------------------------------------------------------------------------------------------------------------------------------------
 */
typedef struct {
        float * o;       // origin
        float * d;       // direction
} Ray;

int normalizeVector3D(float * e)
{
	float magnitude = sqrtf(e[0] * e[0] + e[1] * e[1] + e[2] * e[2]);

	e[0] /= magnitude;
	e[1] /= magnitude;
	e[2] /= magnitude;

	return 0;
}
int normalizeVector4D(float * e)
{
	float magnitude = sqrtf(e[0] * e[0] + e[1] * e[1] + e[2] * e[2] + e[3]*e[3]);

	e[0] /= magnitude;
	e[1] /= magnitude;
	e[2] /= magnitude;
	e[3] /= magnitude;

	return 0;
}

int matrixVector4D(float * matrix, float * v, float * vo)
{
        float x = matrix[0]  * v[0] + matrix[1]  * v[1] + matrix[2]  * v[2] + matrix[3]  * v[3];
        float y = matrix[4]  * v[0] + matrix[5]  * v[1] + matrix[6]  * v[2] + matrix[7]  * v[3];
        float z = matrix[8]  * v[0] + matrix[9]  * v[1] + matrix[10] * v[2] + matrix[11] * v[3];
        float k = matrix[12] * v[0] + matrix[13] * v[1] + matrix[14] * v[2] + matrix[15] * v[3];

	vo[0]=x;
	vo[1]=y;
	vo[2]=z;
	vo[3]=k;

	return 0;
} 

float dot3(float * v1, float * v2)
{
        return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2];
}

float dot4(float * v1, float * v2)
{
        return v1[0] * v2[0] + v1[1] * v2[1] + v1[2] * v2[2] + v1[3] * v2[3];
}

void mul3x4_3(float * M,float * v, float * result)
{
	float x = dot3(v, &M[0]);
    	float y = dot3(v, &M[4]);
    	float z = dot3(v, &M[8]);

    	result[0] = x;
	result[1] = y;
	result[2] = z;   
}

void mul3x4_4(float * M, float * v, float * result)
{
    	float x = dot4(v, &M[0]);
    	float y = dot4(v, &M[4]);
    	float z = dot4(v, &M[8]);

        result[0] = x;
        result[1] = y;
        result[2] = z;
	result[3] = 1.0;   
}


float saturatef(float x){

	if (x < 0.0) 
		return 0.0;
	if (x > 1.0) 
		return 1.0;
	if ((x >= 0.0) && (x <= 1.0 )) 
		return x;
	else 
		return 0.0;
}
uint rgbaFloatToInt(float * rgba)
{
    rgba[0] = saturatef(rgba[0]);   // clamp to [0.0, 1.0]
    rgba[1] = saturatef(rgba[1]);
    rgba[2] = saturatef(rgba[2]);
    rgba[3] = saturatef(rgba[3]);
    return ((uint)(rgba[3]*255)<<24) | ((uint)(rgba[2]*255)<<16) | ((uint)(rgba[1]*255)<<8) | (uint)(rgba[0]*255);
}

int intersecBox(Ray r, float * boxmin, float * boxmax, float *tnear, float *tfar)
{
    // compute intersection of ray with all six bbox planes
	float * invR = (float*) malloc(3*sizeof(float));
	float * tbot = (float*) malloc(3*sizeof(float));
	float * ttop = (float*) malloc(3*sizeof(float));
	float * tmin = (float*) malloc(3*sizeof(float));
	float * tmax = (float*) malloc(3*sizeof(float));
	invR[0] = 1.0 / r.d[0];
	invR[1] = 1.0 / r.d[1];
	invR[2] = 1.0 / r.d[2];
	tbot[0] = invR[0] * (boxmin[0] - r.o[0]);
	tbot[1] = invR[1] * (boxmin[1] - r.o[1]);
	tbot[2] = invR[2] * (boxmin[2] - r.o[2]);
	ttop[0] = invR[0] * (boxmax[0] - r.o[0]);
	ttop[1] = invR[1] * (boxmax[1] - r.o[1]);
	ttop[2] = invR[2] * (boxmax[2] - r.o[2]);

	// re-order intersections to find smallest and largest on each axis
	tmin[0] = fminf(ttop[0],tbot[0]);
	tmin[1] = fminf(ttop[1],tbot[1]);
	tmin[2] = fminf(ttop[2],tbot[2]);
	tmax[0] = fmaxf(ttop[0],tbot[0]);
	tmax[1] = fmaxf(ttop[1],tbot[1]);
	tmax[2] = fmaxf(ttop[2],tbot[2]);

    	// find the largest tmin and the smallest tmax
    	float largest_tmin = fmaxf(fmaxf(tmin[0], tmin[1]), fmaxf(tmin[0], tmin[2]));
    	float smallest_tmax = fminf(fminf(tmax[0], tmax[1]), fminf(tmax[0], tmax[2]));

        *tnear = largest_tmin;
        *tfar = smallest_tmax;

        return smallest_tmax > largest_tmin;
}


void RayCasting(uint *d_output, uint imageW, uint imageH, float density, float brightness, float transferOffset, float transferScale)
{

	int		maxSteps 	= 500;
	int		i		= 0;
	float 		tstep 		= 0.01f;
	float 		opacityThreshold= 0.95f;
	float * 	boxMin 		= (float*) malloc(3*sizeof(float));
	float * 	boxMax 		= (float*) malloc(3*sizeof(float));
	float * 	aux		= (float*) malloc(4*sizeof(float));
	float * 	sum		= (float*) malloc(4*sizeof(float));
	float * 	pos		= (float*) malloc(4*sizeof(float));
	float * 	step		= (float*) malloc(4*sizeof(float));
	float * 	col		= (float*) malloc(4*sizeof(float));
	Ray 		eyeRay;
			eyeRay.o	= (float*) malloc(3*sizeof(float));	
			eyeRay.d	= (float*) malloc(3*sizeof(float));	

	boxMin[0] = -1.0;	boxMin[1] = -1.0;	boxMin[2] = -1.0;
	boxMax[0] = 1.0;	boxMax[1] = 1.0;	boxMax[2] = 1.0;

	uint x = 0;
	uint y = 0;

	for(x=0;x<imageW;x++)
	{
		float u = (x / (float) imageW)*2.0f-1.0f;
		for(y=0;y< imageH;y++)
		{
    			float v = (y / (float) imageH)*2.0f-1.0f;			
			aux[0] = 0.0; aux[1]=0.0; aux[2] = 0.0; aux[3]=1.0;
			mul3x4_4(invViewMatrix,aux,aux);
			eyeRay.o[0] = aux[0]; eyeRay.o[1] = aux[1]; eyeRay.o[2]=aux[2];

			eyeRay.d[0] = u; eyeRay.d[1]=v; eyeRay.d[2] = -2.0;
			normalizeVector3D(eyeRay.d);
			mul3x4_3(invViewMatrix,eyeRay.d,eyeRay.d);
		
			// find intersection with box
			float tnear, tfar;
		        int hit = 0;
			hit = intersecBox(eyeRay, boxMin, boxMax, &tnear, &tfar);
			if (!hit) return;
		        if (tnear < 0.0f) tnear = 0.0f;     // clamp to near plane
			    // march along ray from front to back, accumulating color

			sum[0]=0.0;sum[1]=0.0;sum[2]=0.0;sum[3]=0.0;
		    	float t = tnear;
			pos[0] = eyeRay.o[0] + eyeRay.d[0]*tnear;
			pos[1] = eyeRay.o[1] + eyeRay.d[1]*tnear;
			pos[2] = eyeRay.o[2] + eyeRay.d[2]*tnear;
			step[0]= eyeRay.d[0]*tstep;
			step[1]= eyeRay.d[1]*tstep;
			step[2]= eyeRay.d[2]*tstep;

		    	for(i=0; i<maxSteps; i++) {
			fprintf(stdout,"tocado\n");
			// read from 3D texture
			// remap position to [0, 1] coordinates
			float sample = volume.data[(int)(pos[0]*0.5+0.5)][(int)(pos[1]*0.5+0.5)][(int)(pos[2]*0.5+0.5)];
			//sample *= 64.0f;    // scale for 10-bit data

			// lookup in transfer function texture
			col[0]=pos[0];
			col[1]=pos[1];
			col[2]=pos[2];
			col[3]=0.0;
			normalizeVector4D(col);
			col[3] *= density;

			// "under" operator for back-to-front blending
			//sum = lerp(sum, col, col.w);

			// pre-multiply alpha
			col[0] *= 255;//col.w;
			col[1] *= 255;//col.w;
			col[2] *= 255;//col.w;
			// "over" operator for front-to-back blending
			sum[0] = sum[0] + col[0]*(1.0-sum[3]);
			sum[1] = sum[1] + col[1]*(1.0-sum[3]);
			sum[2] = sum[2] + col[2]*(1.0-sum[3]);
			sum[3] = sum[3] + col[3]*(1.0-sum[3]);
			
			// exit early if opaque
			if (sum[3] > opacityThreshold)
			    break;

			t += tstep;
			if (t > tfar) break;

			pos[0] += step[0];
			pos[1] += step[1];
			pos[2] += step[2];
		    }
		    sum[0] *= brightness;
		    sum[1] *= brightness;
		    sum[2] *= brightness;
		    sum[3] *= brightness;

		    // write output color
		    d_output[y*imageW + x] = rgbaFloatToInt(sum);


		}
	}

	free(boxMin);
	free(boxMax);
	free(aux);
	free(pos);
	free(sum);
	free(step);

}
/*
 *------------------------------------------------------------------------------------------------------------------------------------------------------------------
 *------------------------------------------------------------------------------------------------------------------------------------------------------------------
 */

void initPixelBuffer();

// render image using CUDA
void render()
{

	d_output = realloc(d_output,width*height*sizeof(uint));
	RayCasting(d_output, width,height, density, brightness, transferOffset, transferScale);
	fprintf(stdout,"Renderizar\n");
}

void display()
{

    // use OpenGL to build view matrix
    GLfloat modelView[16];
    glMatrixMode(GL_MODELVIEW);
    glPushMatrix();
        glLoadIdentity();
        glRotatef(-viewRotation[0], 1.0, 0.0, 0.0);
        glRotatef(-viewRotation[1], 0.0, 1.0, 0.0);
        glTranslatef(-viewTranslation[0], -viewTranslation[1], -viewTranslation[2]);
    glGetFloatv(GL_MODELVIEW_MATRIX, modelView);
    glPopMatrix();

    invViewMatrix[0] = modelView[0]; invViewMatrix[1] = modelView[4]; invViewMatrix[2] = modelView[8]; invViewMatrix[3] = modelView[12];
    invViewMatrix[4] = modelView[1]; invViewMatrix[5] = modelView[5]; invViewMatrix[6] = modelView[9]; invViewMatrix[7] = modelView[13];
    invViewMatrix[8] = modelView[2]; invViewMatrix[9] = modelView[6]; invViewMatrix[10] = modelView[10]; invViewMatrix[11] = modelView[14];

    render();

    // display results
    glClear(GL_COLOR_BUFFER_BIT);

    // draw image from PBO
    glDisable(GL_DEPTH_TEST);

    glPixelStorei(GL_UNPACK_ALIGNMENT, 1);
    // draw using texture

    // copy from pbo to texture
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
    glBindTexture(GL_TEXTURE_2D, tex);
    glTexSubImage2D(GL_TEXTURE_2D, 0, 0, 0, width, height, GL_RGBA, GL_UNSIGNED_BYTE, 0);
    glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);

    // draw textured quad
    glEnable(GL_TEXTURE_2D);
    glBegin(GL_QUADS);
    glTexCoord2f(0, 0); glVertex2f(0, 0);
    glTexCoord2f(1, 0); glVertex2f(1, 0);
    glTexCoord2f(1, 1); glVertex2f(1, 1);
    glTexCoord2f(0, 1); glVertex2f(0, 1);
    glEnd();

    glDisable(GL_TEXTURE_2D);
    glBindTexture(GL_TEXTURE_2D, 0);

    glutSwapBuffers();
    glutReportErrors();

}

void idle()
{
    glutPostRedisplay();
}


void keyboard(unsigned char key, int x, int y)
{
    switch(key) {
        case 'f':
            linearFiltering = !linearFiltering;
            break;
        case '+':
            density += 0.01f;
            break;
        case '-':
            density -= 0.01f;
            break;

        case ']':
            brightness += 0.1f;
            break;
        case '[':
            brightness -= 0.1f;
            break;

        case ';':
            transferOffset += 0.01f;
            break;
        case '\'':
            transferOffset -= 0.01f;
            break;

        case '.':
            transferScale += 0.01f;
            break;
        case ',':
            transferScale -= 0.01f;
            break;

        default:
            break;
    }
    glutPostRedisplay();
}


int ox, oy;
int buttonState = 0;

void mouse(int button, int state, int x, int y)
{
    if (state == GLUT_DOWN)
        buttonState  |= 1<<button;
    else if (state == GLUT_UP)
        buttonState = 0;

    ox = x; oy = y;
    glutPostRedisplay();
}


void motion(int x, int y)
{
    float dx, dy;
    dx = (float)(x - ox);
    dy = (float)(y - oy);

    if (buttonState == 4) {
        // right = zoom
        viewTranslation[2] += dy / 100.0f;
    }
    else if (buttonState == 2) {
        // middle = translate
        viewTranslation[0] += dx / 100.0f;
        viewTranslation[1] -= dy / 100.0f;
    }
    else if (buttonState == 1) {
        // left = rotate
        viewRotation[0] += dy / 5.0f;
        viewRotation[1] += dx / 5.0f;
    }

    ox = x; oy = y;
    glutPostRedisplay();
}

int iDivUp(int a, int b){
    return (a % b != 0) ? (a / b + 1) : (a / b);
}

void reshape(int w, int h)
{
    width = w; height = h;
    initPixelBuffer();

    glViewport(0, 0, w, h);

    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();

    glMatrixMode(GL_PROJECTION);
    glLoadIdentity();
    glOrtho(0.0, 1.0, 0.0, 1.0, 0.0, 1.0);
}

void cleanup()
{
    if (pbo) {
        glDeleteBuffersARB(1, &pbo);
        glDeleteTextures(1, &tex);
    }

}

void initGL(int *argc, char **argv)
{
    // initialize GLUT callback functions
    glutInit(argc, argv);
    glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE);
    glutInitWindowSize(width, height);
    glutCreateWindow("CUDA volume rendering");

    glewInit();
}

void initPixelBuffer()
{
    if (pbo) {

        // delete old buffer
        glDeleteBuffersARB(1, &pbo);
        glDeleteTextures(1, &tex);
    }

    // create pixel buffer object for display
    glGenBuffersARB(1, &pbo);
        glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, pbo);
        glBufferDataARB(GL_PIXEL_UNPACK_BUFFER_ARB, width*height*sizeof(GLubyte)*4, 0, GL_STREAM_DRAW_ARB);
        glBindBufferARB(GL_PIXEL_UNPACK_BUFFER_ARB, 0);

    // create texture for display
    glGenTextures(1, &tex);
    glBindTexture(GL_TEXTURE_2D, tex);
    glTexImage2D(GL_TEXTURE_2D, 0, GL_RGBA8, width, height, 0, GL_RGBA, GL_UNSIGNED_BYTE, NULL);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
    glBindTexture(GL_TEXTURE_2D, 0);
}



int main(int argc, char ** argv)
{
	viewRotation    	= (float*) malloc(3*sizeof(float));
	viewRotation[0] = 0.0;	viewRotation[1] = 0.0;	viewRotation[2] = -4.0;
	viewTranslation 	= (float*) malloc(3*sizeof(float));

	initGL( &argc, argv );

	float pos[3] = {0.0, 0.0, 0.0};
	readHDF5(argv[1], "value", &volume);
	setPositionVolume(&volume, pos);

	glutDisplayFunc(display);
        glutKeyboardFunc(keyboard);
        glutMouseFunc(mouse);
        glutMotionFunc(motion);
        glutReshapeFunc(reshape);
        glutIdleFunc(idle);

        initPixelBuffer();

        atexit(cleanup);

        glutMainLoop();

	free(viewRotation);
	free(viewTranslation);

	return 0;
}
